Spin diffusion in trapped clouds of cold atoms with resonant interactions 
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We show that puzzling recent experimental results on spin diffusion in a strongly interacting 
atomic gas may be understood in terms of the predicted spin diffusion coefficient for a generic 
strongly interacting system. Three important features play a central role: a) Fick's law for diffusion 
must be modified to allow for the trapping potential, b) the diffusion coefficient is inhomogeneous, 
due to the density variations in the cloud and c) the diffusion approximation fails in the outer parts 
of the cloud, where the mean free path is long. 



Diffusion in the presence of an external potential is 
an important problem in diverse fields, ranging from 
astrophysics, to condensed matter physics, to biology. 
New vistas for understanding diffusion of spin have been 
opened up by experiments using resonantly interacting 
atomic gases [TJ |2] ■ These experiments are the analog for 
spin phenomena of earlier groundbreaking experiments 
that established that atomic gases may form a perfect 
fluid with a shear viscosity having the least possible value 
consistent with quantum mechanics [3J H] . 

In the spin transport experiments, a cloud of atoms 
consisting of two hyperfine states of the same atom, 
which we refer to as f and J,, was studied. Atoms in 
one state were displaced with respect to those in the 
other state and the subsequent dynamics was investi- 
gated [U [2]- When the population of one hyperfine 
species is much larger than the other, the diffusive mo- 
tion is well described by collisional relaxation [H [5] . For 
equal populations of the f & n d -J. atoms, previous stud- 
ies have focussed on the initial bouncing motion of the 
clouds [fH [7] . Here, we analyze the long time scale dy- 
namics and show that, because of the trap potential V(r), 
Fick's law must be modified. We demonstrate that this 
effect, combined with the fact that the spin diffusion co- 
efficient is inhomogeneous leads to predictions for the de- 
cay rate that are more than f order of magnitude larger 
than the experimentally measured one in the classical 
regime. The resolution of this puzzle is shown to be the 
failure of the diffusion approximation in the outer regions 
of the cloud. Our analysis accounts for the experimen- 
tal results in Ref. pQ using the spin diffusion coefficient 
predicted for a resonantly interacting system. There is a 
rich variety of regimes for spin relaxation, depending on 
the trap anisotropy and the density of atoms. 
Basic formalism In a trap, the magnetization density 
M(r) = n^(r) — n^(r), where «i(r) is the density of 
species i, is not constant in equilibrium. For instance, 
M(r) oc e _v ( r )/ T for high temperatures T. (We use units 
in which ks — 1.) Rather the quantity that is constant 
is the chemical potential difference ^ — /i^ ~ 2M/x, 
where \ = 2dM/d(fif — m) is the spin susceptibility. 
Thus diffusion is driven by spatial variations of the chem- 



ical potentials, and phenomenologically, the spin current 
density jm is therefore given by the modified Fick's law 

hi = -D X V W , (I) 

where D is the spin diffusion coefficient. Equation ([I]) 
reduces to the usual expression jjv/ = — DVM when V 
is constant. We concentrate on the case of tempera- 
tures high enough that the gas may be treated using 
the Maxwell-Boltzmann distribution. In atomic gases, 
the dominant relaxation process is two-body scattering 
and, consequently, the diffusion coefficient, which is pro- 
portional to the mean free path of a particle, therefore 
varies inversely with the density n cx e~ v / T [lj [8] and 
X = n/T, where n = rif + n^. 

To determine the diffusive modes, we write M(r,t) = 
e~ rt M(r). Insertion of Eq. |T]) in the equation of conti- 
nuity, d t M + V • j*m = gives 

A)V 2 P + re" v/T P = 0, (2) 

where Do is the diffusion coefficient at the center of the 
trap (V = 0), and P(r) = M(r)/n(r) is the local frac- 
tional polarization. Equation ^ describes diffusion in 
the presence of an external potential and is often referred 
to as the Smoluchowski equation [9 . In regions where 
V S> T, P satisfies the Laplace equation, and therefore 
the component of P proportional to the spherical har- 
monic Yi m must vary as r~ l in three dimensions since 
the solution varying as r l is forbidden by the condition 
that, by definition, |P| < 1. Thus both P and VP van- 
ish as r — ¥ oo. Equation ^ is therefore analogous to 
the Schrodinger equation for a potential oc e~ v / T and 
determining the eigenvalue T is equivalent to finding the 
strength of the potential that will produce a zero energy 
bound state. Equation ^ may be derived from the vari- 
ational principle <5r var — 0, where 

var ~ JtPre-v/TP 2 ' 1 ' 

For the lowest mode with a particular symmetry, r var 
provides an upper bound on the lowest eigenvalue. 
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FIG. 1: (color on-line) (a) The polarization (solid line) and 
the spin current density (dashed line) for the ID case, (b) 
P(r) for the spherical case with P(r) = P(r) cos 6. 



Simple examples We first solve |2| for a one-dimensional 
(ID) harmonic potential, V — muj 2 z 2 /2, where w z is the 
trap frequency. For \z\ — > oo, P varies as A + Bz, where 
A and B are constants, but because \P\ < 1, B = 0. 
A numerical solution of ^ for the lowest mode that is 
odd in z for the boundary condition P{z) — > constant for 
z — > oo yields for the damping rate the result 



r u 



2.684- 



(4) 



where l\ — 2kBT/mu>f. Plots of the polarization and 
the associated spin current density are given in Fig. [T] 
(a). Since pi) is linear, the normalization of P and 
]m hi Figs, fflp] is arbitrary. The variational function 
tanh(z/0.7842/ z ) gives for r var the value 2.687D /l 2 , 
which is within ~ 0.1% of the exact result. 

We now consider the spherically symmetric case V — 
mtu 2 r 2 /2. The simplest solution rotationally invari- 
ant about the z-axis and odd in z has the form P = 
Y 10 (9)u(r)/r with Y 10 (6) cx cos0. In Fig. [l](b), we plot 
a numerical solution to ([2]). Requiring the solution to 
vanish as r — > oo yields the damping rate 



T = 12.10 



Do 
I 2 ' 



(5) 



Figure [2] (a) shows contour plots of P and the spin cur- 
rent density, which resembles that for a dipole. The vari- 
ational function z/[l + (r/d) 3 } has the correct asymp- 
totic behavior for both r — ¥ and r — > oo, and it yields 
r w = 12.12£>oA 2 for d = 0.886/. 

Anisotropic traps The spin diffusion experiments [1] are 
performed in a prolate trap of the form V(r) = m(w^p 2 + 
oj 2 z 2 )/2 = V± + V z , where p = (x, y), with cjj_ > u». We 
now solve the diffusion equation ^ for a general aspect 
ratio A = aj z /u>± using the variational function 



P(p,z) = z/(l + R 3 ), 



(6) 



with R 2 = p 2 / d 2 L + z 2 / d 2 which obeys the correct bound- 
ary conditions for r —¥ and r — ¥ oo. The variational 
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FIG. 2: (color online) Contour plots of the polarization and 
spin current density (arrows) in the pz-plane. The red dashed 
contour shows where the density has fallen to 0.1 of the central 
value, (a) The spherical case, (b) The prolate case with 
A = uj z /uj±_ = 1/5. (c) The oblate case with A = 5. 
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FIG. 3: (color online) The damping rate as a function of the 
aspect ratio obtained with the variational function Q (solid 
line). The result (JsJ for a spherical t rap is plotted as a cross, 
and the limits A — > oo ^ and A — > ( 11 \ as dashed lines. The 
variational length scales d±_ and d z are shown in the inset. 



parameters d± and d z determine the fall-off of the po- 
larization in the transverse and axial directions in units 
of l± and l z , respectively. The resulting damping rate 
is plotted versus A in Fig. [3] We see that the varia- 
tional function reproduces very accurately the result for 
the spherical case A = 1. The damping is a decreasing 
function of A, since the transverse confinement imposes 
a gradient in the polarization as is illustrated in Fig. [2] 
The inset demonstrates that for prolate traps, the length 
scale of the transverse variations d± becomes longer than 
l± while the scale of axial variations d z becomes shorter 
than l z ; the opposite holds for oblate traps. Figure [2] 
(b) illustrates this important point further for the case 
A = 1/5: the polarization distribution is considerably 
less prolate than the density distribution, and the current 
density is significant even in regions where the density is 
low. Note that, for the prolate and spherical cases, the 
current has large transverse as well as axial components. 

For A — ► oo, we see from Fig. [3] that the damping rate 
approaches the ID result Q. This reflects that the spin 
motion becomes ID with the current essentially in the 
axial direction from the maximum to the minimum of 
the polarization, as is clearly seen in Fig. [2] (c). 

Born-Oppenheimer approximation Using the dimen- 
sionlcss variables p — p/l± and z = z/l z , we see that 
the diffusion equation ^ for an anisotropic, harmonic 
trapping potential is equivalent to a threshold problem 
in quantum mechanics with an isotropic 3D Gaussian 
potential, but where the mass for motion in the trans- 
verse directions is a factor A 2 = uj 2 /u>\ smaller than that 
for axial motion. For the case of a very prolate trap 
with A < 1, we can therefore solve the diffusion equa- 



tion using the Born-Oppenheimer approximation. Writ- 
ing P(r) = tp(z)(j)(p, z), we first find the lowest eigenstate 
for the "light" particle by solving 



dl + dl - A{z)e~~P ct>(p,z) = E(z)<f>(p,z) (7) 



with A{z) = Ao exp(—z 2 ) and Aq 
determined by solving the equation 



Tl 2 JD ; (f>(z) is 



dl + E(z) 



(z) = 0. 



(8) 



The damping T is determined from the value of Aq re- 
quired for Eq. ^ to have a zero-energy bound state that 
is odd under reflection, <p(— z) = —<fr(z). To solve Q, we 
observe that for ujj_ uj z we expect T <C DqIJ 2 , which 
corresponds to Aq <C 1. The transverse problem then re- 
duces to finding the energy of the lowest bound state in a 
shallow 2D Gaussian potential V(p) = —A(z)exp(—p 2 ). 
Such a state always exists and for Vq <C 1 its energy is 
E = -Kexp(47r/F ) with V = J d 2 pV(p) QUE]. The 
prefactor is K — 2 exp(— 27 + D\) to lowest order in A(z) 
where 7 ~ 0.5772 is Euler's constant and [T2] 



D 1 = / d z pxd 2 p 2 



s V(pi)V(p 2 



ln(p 2 2 /2) (9) 



V 2 
1 



with pi 2 = pi — f>2- The integrals are straightforward to 
perform, and we find 



E(z) = -2e 



-3 7e -4/A(z)_ 



(10) 



For Aq -c 1, we can expand the eigenvalue as E(z) ~ 
— 2exp[— 37 — 4(1 + z 2 )/Aq] and when this is inserted in 
([8]), we recover the ID diffusion problem in a Gaussian 
trap. Using our ID result Q, the threshold condition for 
a bound state odd in z in a very prolate trap becomes 



Te 



-4D /r/i _ 2e 37p 



ID, 



(11) 



which is an implicit equation for the damping rate. For 
T <C 1, the solution is 



4 Dp 
\n(P z /7.6l 2 ± ) IY 



(12) 



Thus, the assumption Aq <C 1 for l± <C l z is consistent. 
However, the numerical factors indicate that the asymp- 



totic expression ( 12 ) is a good approximation only for 



extremely prolate traps. Figure [3] shows that the varia- 
tional function ^ accurately recovers the prolate limit 
of the damping rate given by ( |11| . 

Failure of the diffusion approximation For A w 0.1 
the calculations above predict a damping rate of ap- 
proximately 200-DoAz, while experimentally the results 
of Ref. 4 with the expression for D from Ref. [5j give 
r m IOD0/I 2 . We now demonstrate that the discrepancy 
is due to the failure of the diffusion approximation in the 
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outer parts of the cloud, where the density is low and 
conditions are collisionless. An approximate expression 
for the distance r$ from the z axis at which the diffu- 
sion approximation fails may be obtained by arguing that 
this occurs when a particle has a probability of l/e" of 
not suffering a collision when it comes in from infinity. 
Assumingr ~> l±, this gives [13] 



(13) 



along the z-direction with a damping rate given by ( 14 1 




Equation ( 13 1 is correct to logarithmic accuracy. The 



exact value of the parameter a ~ 0(1) may be deter- 
mined by solving the kinetic equation in the vicinity of 
the boundary. However, provided tq 3> l±_ the uncer- 
tainty in a has little effect on the value of r$. 

At p — tq it is necessary to impose a boundary con- 
dition. The flux of atoms in the p-direction for p > r$ 
is small since, in the absence of collisions, atoms mov- 
ing to larger values of p will be reflected by the trapping 
potential thereby strongly reducing the net flux. Conse- 
quently, we impose the condition that the current in the 
p-direction vanish at p = ro, or dP/dp\ p=ro = 0. When 
the distance to the boundary ro is not much larger than 
the typical length scale I for the spin diffusion modes, 
it influences the damping rate. In a prolate trap for 
which r > £j_> the polarization reaches its asymptotic 
for z — > oo for \z\ < l z and consequently the failure of the 
diffusion approximation for the motion in the axial direc- 
tion, which will occur only for distances large compared 
with l z , does not affect the results. With the boundary 
condition dP/dp\ p=ra — 0, the variational principle de- 
rived earlier applies except that the region of integration 
is limited to p < tq. For ro not too much larger than 
we argue that a good approximation for a trial function 
is simply a function of z, which gives a damping rate 



IT 



(14) 



With increasing r , there will be a cross-over between 
the ID spin currents with a damping given by ( fl4| ) and 
the fully 3D hydrodynamic spin currents which have both 
transverse and axial components as shown in Fig. [I] (a)- 
(b) with a damping scaling as T ~ D /l\. We expect 
the cross-over between the two hydrodynamic solutions 
to occur when ro/J 2 r 1D ~ D /l\ which gives r ~ l z . 

Comparison with experiment We finally compare our 
results with the experiments in Ref. [T]. For long times, 
the spin dynamics is determined by the lowest diffu- 
sive mode and the observed decay time is related to 
the damping rate by the relation r = 1/T. Using typ- 
ical experimental numbers reported in Ref. [1], we ob- 
tain 3 < < 6. This means that the spin dynam- 
ics is diffusive in a major part of the cloud, and since 
l z /l±_ ~ 10 > ro/l± we expect the motion to be mainly 



In Ref. pP, the decay time is written in terms of a spin 
drag coefficient as r — T s d/oj 2 and we find 



nr. 



sd 



Et 



I 2 



Dn 



h T 



I T m l 2 z mD Tp 



I 2 



(15) 



where we have used (14 1 and D ~ l.l(T/T F ) 3 / 2 h/m 



for a strongly interacting Fermi gas in the classical 
regime (TJ [5] . This agrees with the measured high tem- 
perature experimental result, T sc i = OA&Epfrr 1 \jTp /T 
when 7gj_/Zj_ ~ 4.4 which is consistent with the estimate 
above. In addition to reproducing the magnitude and 
temperature dependence of the damping, our result also 
explains the observation that T s d is independent of the 
axial trapping frequency uj z . 

In summary, we have shown that a quantitative ac- 
count of the measured damping rates of diffusive modes 
can be given if three novel features of spin diffusion in a 
trap are taken into account: the failure of Fick's law, the 
inhomogeneity of the diffusion coefficient, and the failure 
of the diffusion approximation in the outermost regions 
of the cloud. The work may be extended in a number 
of directions: to gases with unequal numbers of the two 
species, and to degenerate gases, including ones with a 
condensate of paired fermions. Calculations of damping 
rates may also be refined by improving variational func- 
tions and by obtaining a more quantitative understand- 
ing of the boundary between diffusive and collisionless 
behavior on the basis of the Boltzmann equation. 
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